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Orbital-free (OF) methods promise significant speed-up of computations based 
on density functional theory (DFT). In this field, the development of accurate 
kinetic-energy density functionals remains an open question. In this chapter 
we review the shell-correction method (SCM, commonly known as Strutinsky's 
averaging method) applied originally in nuclear physics and its more recent for- 
mulation in the context of DFT [Yannouleas and Landman, Phys. Rev. B 48, 
8376 (1993)]. We demonstrate the DFT-SCM method through its earlier applica- 
tions to condensed-matter finite systems, including metal clusters, fuUerenes, and 
metal nanowires. The DFT-SCM incorporates quantum mechanical interference 
effects and thus offers an improvement compared to the use of Thomas- Fermi-type 
kinetic energy density functionals in OF-DFT. 

1.1. Introduction 
1.1.1. Preamble 

Often theoretical methods (in particular computational techniques) are developed in 
response to emerging scientific challenges in specific fields. The development of the 
shell correction method (SCM) by Strutinsky in the late 1960's [1] was motivated 
by the observation of large nonuniformities (oscillatory behavior) exhibited by a 
number of nuclear properties as a function of the nuclear size. These properties 
included: total nuclear masses, nuclear deformation energies, and large distortions 
and fission barriers. While it •was understood already that the total energy of nuclei 
can be decomposed into an oscillatory part and one that sho"ws a slo"w "smooth" 
variation as a function of size, Strutinsky's seminal contribution was to calculate the 
two parts from different nuclear models: the former from the nuclear shell model 
and the latter from the liquid drop model. In particular, the calculation of the 
oscillatory part was enabled by employing an averaging method that smeared the 
single particle spectrum associated with a nuclear model potential. It is recognized 
that the Strutinsky procedure provides a method which "reproduces microscopic 
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results in an optimal 'way using phenomenological models" [2]; in the Appendix 
we describe an adaptation of the Strutinsky phenomenological procedure to metal 
clusters; we term this procedure as the semi-empirical (SE)-SCM. 

In this chapter, we focus on our development in the early 1990's of the micro- 
scopic density- functional-theory (DFT)-SCM approach [3], where we have shown 
that the total energy of a condensed-matter finite system can be identified with 
the Harris functional [see Eq. (|1.16p ]. with the shell correction [Eq. (|1.23|) ] being 
expressed through both the kinetic energy, Tgh, of this functional [Eq. (|1.19l) ] and 
the kinetic energy of an extended-Thomas-Fermi (ETF) functional expanded to 
fourth-order density gradients [see Tetf in Eq. (|1.22p ]. It is important to note that 
in our procedure an optimized input density is used in the Harris functional. This 
optimization can be achieved through a variational procedure [using an orbital-free 
(OF) energy functional, e.g., the ETF functional with 4th-order gradients] with 
a parametrized trial density profile [see Eq. (|1.25p ]. or through the use [see Eq. 
(|1.24[) ] of the variational principle applied to an orbital-free energy functional. (For 
literature regarding orbital-free kinetic-energy functional, see, e.g., Refs. [4-9].) A 
similar optimization of an OF/4th-order-ETF density has been shown to be con- 
sistent with the Strutinsky averaging approach [10]. Such 4th-ordcr optimization 
of the input density renders rather ambiguous any direct (term-by-term) compari- 
son between the method proposed by us and subsequent treatments, which extend 
the DFT-SCM to include higher-order shell-correction terms without input-density 
optimization [11, 12] (see also Ref. [13]). Indeed, the input-density optimization 
(in particular with the use of 4th-order gradients) minimizes contributions from 
higher-order shell corrections. 

In light of certain existing similarities between the physics of nuclei and clusters 
(despite the large disparity in spatial and energy scales and the different origins of 
inter-particle interactions in these systems), in particular the finding of electronic 
shell effects in clusters [14-18], it was natural to use the jellium model in the early 
applications of the DFT-SCM to clusters. However, as noted [3] already early on, 
"the very good agreement between our results and those obtained via Kohn-Sham 
self-consistent jellium calculations suggested that it would be worthwhile to extend 
the application of our method to more general electronic structure calculations ex- 
tending beyond the jellium model, where the trial density used for minimization of 
the ETF functional could be taken as a superposition of site densities, as in the 
Harris method." Additionally, generalization of the DFT-SCM method to calcu- 
lations of extended (bulk and surface) systems appeared rather natural. Indeed, 
recent promising applications of DFT-SCM in this spirit have appeared [19, 20]. In 
this case, the term "shell correction effects" is also maintained, although "quantum 
interference effects" could be more appropriate for extended systems. 
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1.1.2. Motivation for finite systems 

One of the principle themes in research on finite systems (e.g., nuclei, atomic and 
molecular clusters, and nano-structured materials) is the search for size-evolutionary 
patterns (SEPs) of properties of such systems and elucidation of the physical prin- 
ciples underlying such patterns [21]. 

Various physical and chemical properties of finite systems exhibit SEPs, includ- 
ing: 

1. Structural characteristics pertaining to atomic arrangements and particle 
morphologies and shapes; 

2. Excitation spectra involving bound-bound transitions, ionization potentials 
(IPs), and electron affinities (EAs); 

3. Collective excitations (electronic and vibrational); 

4. Magnetic properties; 

5. Abundance spectra and stability patterns, and their relation to binding 
and cohesion energetics, and to the pathways and rates of dissociation, 
fragmentation, and fission of charged clusters; 

6. Thermodynamic stability and phase changes; 

7. Chemical reactivity. 

The variations 'with size of certain properties of materials aggregates are com- 
monly found to scale on the average 'with the surface to volume ratio of the cluster, 
i.e., S/fl ^ R^^ ^ N^/^, 'where S, fl, R, and N are the surface area, volume, aver- 
age radius, and number of particles, respectively (the physical origins of such scaling 
may vary for different properties). In general, the behavior of SEPs in finite sys- 
tems in terms of such scaling is non-universal, in the sense that it is non-monotonic 
exhibiting characteristic discontinuities. Nevertheless, in many occasions, it is con- 
venient to analyze the energetics of finite systems in terms of t'wo contributions, 
namely, (i) a term 'which describes the energetics as a function of the system size in 
an average sense (not including shell-closure effects) , referrred to usually as describ- 
ing the "smooth" part of the size dependence, and (ii) an electronic shell-correction 
term. The first term is the one 'which is expected to vary smoothly and be ex- 
pressible as an expansion in S/fl, 'while the second one contains the characteristic 
oscillatory patterns as the size of the finite system is varied. Such a strategy has 
been introduced [1] and often used in studies of nuclei [2, 22], and has been adopted 
recently for investigations of metal clusters [3, 23-38], fullerenes [39], and metal 
nano'wires [40-42]. As a motivating example, we sho'w in Fig. II. li the SEP of the IPs 
of Najv clusters, -which illustrates odd-even oscillations in the observed spectrum, a 
smooth description of the pattern [Fig. Il.lf a)]. and t'wo levels of shell-corrected de- 
scriptions — one assuming spherical symmetry [Fig. ll.ir b)]. and the other allo-wing 
for triaxial shape deformations [Fig. ILlf c)]. The progressive improvement of the 
level of agreement bet'ween the experimental and theoretical patterns is evident. 
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Fig. 1.1. Ionization potentials of Najv clusters. Open squares: Experimental measurements of 
Refs. [15, 17]. The solid line at the top panel (a) represents the smooth contribution to the 
theoretical total IPs. The solid circles in the middle (b) and bottom (c) panels are the total SCM 
IPs. The shapes of sodium clusters have been assumed spherical in the middle panel, while triaxial 
deformations have been considered at the bottom one. 



1.1.3. Plan of the chapter 

The chapter is organized as follows: 

In Sec. II. 2[ the general methodology of shell-correction methods is reviewed, 
and the microscopic DFT-SCM is introduced and presented in detail in Sec. 11.2.21 

Applications of the DFT-SCM to condensed-matter finite systems are presented 
in Sec. 11.31 for three different characteristic nanosy stems, namely, metal clusters 
(Sec. I1.3.T|) . charged fullerenes (Sec. 11.3.21) . and metallic nanowires (Sec. I1.3.3|) . 

In the Appendix, we describe the semiempirical SCM for clusters, which is closer 
to the spirit of Strutinsky's original phenomenological approach for nuclei. There 
we also briefiy present applications of the SE-SCM to triaxial shape deformations 
and fission of metal clusters. 

A summary is given in Sec. 11.41 
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1.2. Methodology and derivation of microscopic DFT-SCM 
1.2.1. Historical review of SCM 

It has long been recognized in nuclear physics that the dependence of ground- 
state properties of nuclei on the number of particles can be vie'wed as the sum of 
t"wo contributions: the first contribution varies smoothly "with the particle number 
(number of protons Np and neutrons N^) and is referred to as the smooth part; 
the second contribution gives a superimposed structure on the smooth curve and 
exhibits an oscillatory behavior, -with extrema at the nuclear magic numbers [22, 43]. 

Nuclear masses have provided a prototype for this behavior [43]. Indeed, the 
main contributions to the experimental nuclear binding energies are smooth func- 
tions of the number of protons and neutrons, and are described by the semi-empirical 
mass formula [44, 45]. The presence of these smooth terms led to the introduction of 
the liquid-drop model (LDM) , according to which the nucleus is viewed as a drop of 
a nonviscous fluid whose total energy is specified by volume, surface, and curvature 
contributions [22, 43, 46]. 

The deviations of the binding energies from the smooth variation implied by the 
LDM have been shown [1, 46] to arise from the shell structure associated with the 
bunching of the discrete single-particle spectra of the nucleons, and are commonly 
referred to as the shell correction. Substantial progress in our understanding of 
the stability of strongly deformed open-shell nuclei and of the dynamics of nuclear 
fission was achieved when Strutinsky proposed [1] a physically motivated efficient 
way of calculating the shell corrections. The method consists of averaging [see the 
Appendix, Eq. (lA.ll) and Eq. (jA.2[) ] the single-particle spectra of phenomenological 
deformed potentials and of subtracting the ensuing average from the total sum of 
single-particle energies. 

While certain analogies, portrayed in experimental data, between properties of 
nuclei and elemental clusters have been recognized, the nuclear-physics approach of 
separating the various quantities as a function of size into a smooth part and a shell 
correction part has only partially been explored in the case of metal clusters. In par- 
ticular, several investigations [47-50] had used the ETF method in conjunction with 
the jellium approximation to determine the average (smooth, in the sense defined 
above) behavior of metal clusters, but had not pursued a method for calculating 
the shell corrections. 

In the absence of a method for appropriately calculating shell-corrections for 
metal clusters in the context of the semiclassical ETF method, it had been presumed 
that the ETF method was most useful for larger clusters, since the shell effects 
diminish with increasing size. Indeed, several studies had been carried out with 
this method addressing the asymptotic behavior of ground-state properties towards 
the behavior of a jellium sphere of infinite size [51, 52]. 

It has been observed [48, 53-56], however, that the single-particle potentials 
resulting from the semiclassical method are very close, even for small cluster sizes. 
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to those obtained via self-consistent solution of the local density functional ap- 
proximation (LDA) using the Kohn-Sham (KS) equations [57]. These semiclassical 
potentials "were used extensively to describe the optical (linear) response of spheri- 
cal metal clusters, for small [53-55], as "well as larger sizes [56] (for an experimental 
revie'w on optical properties, cf. Refs. [58, 59]). The results of this approach are 
consistent "with time-dependent local density functional approximation (TDLDA) 
calculations "which use the KS solutions [60, 61]. 

It is natural to explore the use of these semiclassical potentials, in the spirit of 
Strutinsky's approach, for evaluation of shell corrections in metal clusters of arbi- 
trary size. Belo"w "we describe a microscopic derivation of an SCM in conjunction 
"with the density functional theory [3, 23, 24], and its applications in investiga- 
tions of the properties of metal clusters and fullerenes. Particularly interesting 
and promising is the manner by "which the shell corrections are introduced by us 
at the microscopic level through the kinetic energy term [3, 23, 24], instead of the 
traditional semiempirical Strutinsky averaging procedure of the single-particle spec- 
trum [1, 31]. In particular, our approach leads to an energy functional that corrects 
many shortcomings of the orbital-free DFT, and one that is competitive in numeri- 
cal accuracy and largely advantageous in computational speed compared to the KS 
method. 

1.2.2. DFT-SCM 

Underlying the development of the shell-correction method is the idea of approxi- 
mating the total energy -Etotail-^) of a finite interacting fermion system as 



■where E is the part that varies smoothly as a fimction of system size, and AEsh is an 
oscillatory term. Various implementations of such a separation consist of different 
choices and methods for evaluating the t"wo terms in Eq. (|l.ip . Before discussing 
such methods, "we outline a microscopic derivation of Eq. (jl.ip . 

Motivated by the behavior of the empirical nuclear binding energies, Strutinsky 
conjectured that the self-consistent Hartree-Fock density phf can be decomposed 
into a smooth density p and a fluctuating contribution dp, namely phf = P + 5 p. 
Then, he proceeded to sho"w that, to second-order in 5p, the Hartree-Fock energy 
is equal to the result that the same Hartree-Fock expression yields "when phf is 
replaced by the smooth density p and the Hartree-Fock single-particle energies ef^ 
are replaced by the single-particle energies corresponding to the smooth potential 
constructed "with the smooth density p. Namely, he sho"wed that 



EtoUN) = E{N) + /\E,h{N), 



(1.1) 



EHF^Estr+OiSp"), 

"where the Hartree-Fock electronic energy is given by the expression 



(1.2) 



occ ^ „ 

Ehf ^^^r^ ~2 J 



dvdv'V{v-v')[pHF{v,v)pHF{v'y)~ pHFivyfl (1.3) 
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•with e^^ being the eigenvahies obtained through a self-consistent solution of the 
HF equation, 

-|^V2 + C/ff^.)</),=ef^0„ (1.4) 
•where 



C/i/F(r)0,(r) = j dr'Vir ~ r')[pHFir' ,r')Mr) - PHF{r' ,r)Mr')]- (1-5) 
The Strutinsky approximate energy is written as follows, 

occ ^ „ 

Estr = J2^^-^ / drdr'V(r-r')[p(r,r)p(r',r') -p(r,r')2], (1.6) 
1=1 

■where the index i in Eq. (|1.3p and Eq. (|1.6p runs only over the occupied states (spin 
degeneracy is naturally implied). The single-particle energies correspond to a 
smooth potential U. Namely, they are eigenvalues of a Schrodinger equation, 

— V^ + C/ =£,(p„ (1.7) 



2m 

■where the smooth potential U depends on the smooth density p, i.e., 

C/(r)(p,(r)= /dr'V(r-r')[p(r',r')(p.(r)-p(r',r)^,(r')], (1 



and V is the nuclear t^wo-body interaction potential. 

It should be noted that ■while Eq. (II. 6^ — Eq. (|1.8I) look formally similar to the 
Hartree-Fock equations (|1.3m.5p . their content is different. Specifically, ■while in the 
HF equations, the density phf is self-consistent ■with the ■wavefunction solutions of 
Eq. (|1.4I) , the density p in Eq. (|1.6p — Eq. (11.81) is not self-consistent ■with the ■wave- 
function solutions of the corresponding single-particle Eq. (|1.7p . i.e., p Y^'i=i I'PiP- 
We return to this issue belo^w. 

Since the second term in Eq. (|1.6I) is a smooth quantity, Eq. (|1.2p states that 
all shell corrections are, to first order in 6p, contained in the sum of the single- 
particle energies Y^°^i£i- Consequently, Eq. (|1.6p can be used as a basis for a 
separation of the total energy into smooth and shell-correction terms as in Eq. 
(jl.l|) . Indeed Strutinsky suggested a semi-empirical method of such separation 
through an averaging procedure of the single-particle energies in conjunction 
■with a phenomenological (or semi-empirical) model [the liquid drop model (LDM)] 
for the smooth part (see the appendix). 

Motivated by the above considerations, we have extended them [3, 23, 24] in the 
context of density functional theory for electronic structure calculations. First we 
revie^w pertinent aspects of the DFT theory. In DFT, the total energy is given by 

E[p] = T[p] -f 11 ivH[p(r)] + Vj{r) p(r)| dr + J £.e[p(r)]rfr + Ej, (1.9) 
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■where Vh is the Hartree repulsive potential among the electrons, Vj is the interaction 
potential bet'ween the electrons and ions, £xc is the exchange-correlation functional 
[the corresponding xc potential is given as T4c(r) = S£xcP{t^)/^p{t^)] and T[p] is 
given in terms of a yet unknown functional t[p{r)] as T[p] = J t[p{r)]dr. Ej is the 
interaction energy of the ions. 

In the Kohn-Sham (KS)-DFT theory, the electron density is evaluated from the 
single-particle "wave functions (j)KS,i{^) as 

occ 

PKs(r)=5]|^K5,.(r)|', (1.10) 

1=1 

•where <j)KS,i{'<^) are obtained from a self-consistent solution of the KS equations, 

_9 1 

<t)KS,i{'^) ^ SksAksA'^) (1-11) 



2 m 



•where 

VKs[pKs{r)] = Vh[pks{v)] + 14c[pKs(r)] + V/(r). (1.12) 
The kinetic energy term in Eq. (jl.9D is given by 

occ ^2 

T[pKs\ ^Y.< <^^s,d - :^y^\cbKs,^ >, (1.13) 



2m 

i=l 



■which can also be "written as 

occ 



T[PKS] = ^ei^s.^ - / PKs{r)VKs[pKs(r)]dr. (1.14) 

i—l 



According to the Hohenberg-Kohn theorem, the energy functional p.9p is a 
minimum at the true ground density pgs, ■which in the context of the KS-DFT 
theory corresponds the the density, pks^ obtained from an iterative self-consistent 
solution of Eq. (jl.lip . In other -words, combining Eq. (|1.9p and Eq. (|1.14|) . and 
denoting by "in" and "out" the trial and output densities of an iteration cycle in 
the solution of the KS equation [Eq. ()l.lip ]. one obtains. 



/{ 



j^//[PKs(r)] + ^APkW)\ + V/(r)| pTs(y)dv - 

JpTsi^)VKs[p']isi^)]dr- (1.15) 

Note that the expression on the right involves both p'^g and p^j^s- Self-consistency 
is achieved -when ^p^^'™(r) = P/f 5 (r) — P^s becomes arbitrarily small (i.e., -when 
P'ks converges to pks)- 

On the other hand, it is desirable to introduce approximate energy functionals 
for the calculations of ground-state electronic properties, providing simplified, yet 
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accurate, computational schemes. It is indeed possible to construct such functionals 
[62-66], an example of "which was introduced by J. Harris [62], where self-consistency 
is circumvented and the result is accurate to second order in the difference between 
the trial and the self-consistent KS density (see in particular Eq. (24a) of Ref. [66]; 
the same also holds true for the difference between the trial and the output densities 
of the Harris functional). 

The expression of the Harris functional is obtained from Eq. (ll.lSp by dropping 
the label KS and by replacing everywhere by p'", yielding [note cancellations 
between the third and fourth terms on the right-hand-side of Eq. dl.lSp ]. 

EH.rris[n = Ej + EC* - / [Ivnirir)] + F.c[p"(r)]| p'"(r)dr -|- 

J f,,[p'"(r)]dr. (1.16) 

are the single-particle solutions (non-self-consistent) of Eq. (jl.lip . with 
T^Ks[p"(r) [seeEq. dnil)]. 

As stated above this result is accurate to second order in p^'^ — pks (alternatively 
in p'" — p°"*), thus approximating the self-consistent total energy EksIpks]- 

Obviously the accuracy of the results obtained via Eq. (I1.16P depend on the 
choice of the input density p'". In electronic structure calculations where the cor- 
puscular nature of the ions is included (i.e., all-electron or pseudo-potential calcula- 
tions), a natural choice for p™ consists of a superposition of atomic site densities, as 
suggested originally by Harris. In the case of jellium calculations, we have shown [3] 
that an accurate approximation to the KS-DFT total energy is obtained by using 
the Harris functional with the input density, p'", in Eq. (jl.l6p evaluated from an 
Extended-Thomas-Fermi (ETF)-DFT calculation. 

The ETF-DFT energy functional, Eetf[p], is obtained by replacing the kinetic 
energy term in Eq. (|1.9[) by a kinetic energy density-functional in the spirit of 
the Thomas-Fermi approach [67], but comprising terms up to fourth-order in the 
density gradients [48, 68]. The optimal ETF-DFT total energy is then obtained 
by minimization of Eetf[p\ with respect to the density. In our calculations, we 
use for the trial densities parametrized profiles p(r; {7^}) with {7^} as variational 
parameters (the ETF-DFT optimal density is denoted as petf)- The single-particle 
eigenvalues, in Eq. ()1.16p are obtained then as the solutions to a single- 

particle Hamiltonian, 

Hetf = -t^V^ + Vetf, (1.17) 
Zm 

where Vetf is given by Eq. (|1.12l) with pks{'^) replaced by petf{i^)- These single- 
particle eigenvalues will be denoted by {sj} 

As is well known, the ETF-DFT does not contain shell effects [48-50]. Conse- 
quently, the corresponding density petf can be taken a la Strutinsky as the smooth 
part, p, of the KS density, pks- Accordingly, Eetf is identified with the smooth 
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part E in Eq. (11.11) (in the folio-wing, the "ETF" subscript and " ~ " can be used 
interchangeably). Since, as aforementioned, -Enanis [pbtf] approximates well [i.e., 
to second order in {petf — Pks)] the self-consistent total energy EksIpks], it fol- 
lo^ws from Eq. (|1.1|) . -with i?Harris [pbtf] taken as the expression for -Btotai, that the 
shell-correction, AEsh, is given by 

AEsh - ^HarrisbsTF] " EeTf[PETf] = E,h[p\ " E[p\. (1.18) 

Defining, 



OCC „ 

J PETF{r)VETF{r)dr, (1.19) 



and denoting the total energy E'Harris by Esh, i-C, by identifying 

Esh = i^Harris, (1-20) 

we obtain 

E,h[p\ = {Tsh-f[p\} + E[p\, (1.21) 

■where T[p\ is the ETF kinetic energy, given to fourth-order gradients by the expres- 
sion [68], 



rr-i 1 (Vp)2 j_ 

36 p ^ 27 

4 ^/,-^\2a /a\2 



36 p 270 



— (37r2)-2/3pi/3 



dr, (1.22) 



■which as noted before does not contain shell effects. Therefore, the shell correction 
term in Eq. (|l.ip [or Eq. (|1.18p ] is given by a difference bet^ween kinetic energy 
terms, 

AE,h = Tsh - f[p\. (1.23) 

One should note that the above derivation of the shell correction does not involve 
a Strutinsky averaging procedure of the kinetic energy operator. Rather it is based 
on using ETF quantities as the smooth part for the density, p, and energy, E. Other 
descriptions of the smooth part may result in different shell-correction terms. 

To check the accuracy of this procedure, we have compared results of calculations 
using the functional Esh [Eq. (|1.2ip ] ■with available Kohn-Sham calculations. In 
general, the optimized density from the minimization of the ETF-DFT functional 
can be obtained numerically as a solution of the differential equation 

+ y,,Wp(r)] - M, (1.24) 

■where p, is the chemical potential. As mentioned already, for the jellium DFT-SCM 
calculations, we often use a trial density profile in the ETF-DFT variation which is 
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Fig. 1.2. Total energy per atom of neutral sodium clusters (in units of the absolute value of 
the energy per atom in the bulk, |ecx)| = 2.252 eV). Solid circles; DFT-SCM results (see text for 
details). The solid line is the ETF result (smooth contribution). In both cases, a spherical jellium 
background was used. Open squares: Kohn-Sham DFT results from Ref. [69]. The excellent 
agreement (a discrepancy of only 1%) between the DFT-SCM and the Kohn-Sham DFT approach 
is to be stressed. 



•with ro, a, and 7 as variational parameters that minimize the ETF-DFT functional 
(for other closely related parametrizations, cf. Refs. [49, 50]). 

Figure 11.21 displays results of the present shell correction approach for the total 
energies of neutral sodium clusters. The results of the shell correction method for 
ionization potentials of sodium clusters are displayed in Fig. 11.31 The excellent 
agreement between the oscillating results obtained via our DFT-SCM theory and 
the Kohn-Sham results (cf., e.g., Ref. [69]) is evident. To further illustrate the two 
components (smooth contribution and shell correction) entering into our approach, 
we also display the smooth parts resulting from the ETF method. (In all calcula- 
tions, the Gunnarsson-Lundqvist exchange and correlation energy functionals were 
used; see Refs. [3, 23].) 

1.3. Applications of DFT-SCM 

1.3.1. Metal clusters 

1.3.1.1. Charging of metal clusters 

Investigations of metal clusters based on DFT methods and self-consistent solu- 
tions of the Kohn-Sham equations (employing either a positive jellium background 
or maintaining the discrete ionic cores) have contributed significantly to our un- 
derstanding of these systems [69-72]. However, even for singly negatively charged 
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(1.25) 
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Fig. 1.3. Ionization potentials for sodium clusters. Solid circles: IPs calculated with the DFT- 
SCM (see text for details). The solid line corresponds to the ETF results (smooth contribution). 
In both cases, a spherical jellium background was used. Open squares: Kohn-Sham DFT results 
from Ref. [69]. The excellent agreement (a discrepancy of only 1-2%) between the DFT shell 
correction method and the full Kohn-Sham approach should be noted. 

metal clusters {M^ ) , difBculties may arise due to the failure of the solutions of the 
KS equations to converge, since the eigenvalue of the excess electron may iterate to 
a positive energy [73] . While such difficulties are alleviated for clusters via self- 
interaction corrections (SIC) [74, 75], the treatment of multiply charged clusters 
{M^~ , Z > 1) -would face similar difficulties in the metastability region against 
electronic autodetachment through a Coulombic barrier. In the follo^wing we are 
applying our DFT-SCM approach, described in the previous section, to these sys- 
tems [3, 23, 24] (for the jellium background, we assume spherical symmetry, unless 
otherwise stated; for a discussion of cluster deformations, see Sec. lA.ip . 

1.3.1.2. Electron affinities and borders of stability 

The smooth multiple electron affinities Az prior to shell corrections are defined as 
the difference in the total energies of the clusters 



■where N is the number of atoms, v is the valency and Z is the number of excess 
electrons in the cluster (e.g., first and second affinities correspond to Z — 1 and 
Z ~ 2, respectively). vN is the total charge of the positive background. Applying 
the shell correction in Eq. (I1.23|) . we calculate the full electron afiinity as 



A positive value of the electron affinity indicates stability upon attachment of 
an extra electron. Figure 11.41 displays the smooth, as well as the shell corrected, 
first and second electron affinities for sodium clusters with TV < 100. Note that 



Az = E{vN, vN + Z -1)- E{vN, vN + Z), 



(1.26) 



Af -Az^ AEshivN, vN + Z-l)- AEsh{vN, vN + Z). 



(1.27) 
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Fig. 1.4. Calculated first (Ai) and second {A2) electron affinities of sodium clusters as a function 
of the number of atoms N. Both their smooth part (dashed lines) and the shell-corrected affinities 
(solid circles) are shown. A spherical jellium background was used. 



A2 becomes positive above a certain critical size, implying that the second electron 



in doubly negatively charged sodium clusters -with N < 



(2) 



43 might not be 



stably attached. The shell effects, hcwever, create t-wo islands of stability about the 
magic clusters Nag^ and Nag^ (see Aj'' in Fig. II. 4[) . To predict the critical cluster 



( z) 

size NcT , 'which allo'ws stable attachment of Z excess electrons, we calculated the 
smooth electron affinities of sodium clusters up to = 255 for 1 < Z < 4, and 



display the results in Fig.O We observe that = 205, -while N^^' > 255. 

The similarity of the shapes of the curves in Fig. 11.51 and the regularity of 
distances bet-ween them, suggest that the smooth electron affinities can be fitted by 
a general expression of the form; 



r(4) 



R + S 



[Z-l) 



(1.28) 



R+5 R+S ' 

r,7Vi/3. From our fit, we find 



'where the radius of the positive background is R 
that the constant W corresponds to the bulk -work function. In all cases, we find 
/? = 5/8, 'which suggests a close analogy 'with the classical model of the image 
charge [76, 77]. For the spill-out parameter, we find a -weak size dependence as (5 = 
<5o + 62/ R^- The contribution of S2/R'^, 'which depends on Z, is of importance only 
for smaller sizes and does not affect substantially the critical sizes ('where the curve 



April 21, 2010 0:52 



World Scientific Review Volume - 9.75in x 6. Sin 



yannouleas'of'dft 



14 



Constantine Yannouleas and Uzi Landman 




Fig. 1.5. Calculated smooth electron affinities Az, Z = 1-4, for sodium clusters as a function 
of the number of atoms N (Z is the number of excess electrons). A spherical jellium background 
was used. Inset: The electron drip line for sodium clusters. Clusters stable against spontaneous 
electron emission are located above this line. While for spherical symmetry, as seen from Fig. 11.41 
shell effects influence the border of stability, shell-corrected calculations [25] including deformations 
(see the appendix) yield values close to the drip line (showm in the inset) wrhich was obtained from 
the smooth contributions. 



crosses the zero line), and consequently S2 can be neglected in such estimations. 
Using the values obtained by us for Ai of sodium clusters (namely, W — 2.9 eV 
which is also the value obtained by KS-DFT calculations for an infinite planar 
surface [78], So = 1.16 a.u.; with R ~ VgN^/^, and Ts — 4.00 a. it.), we find for 
the critical sizes when the l.h.s. of Eq. (|1.28p is set equal to zero, Ni? = 44, 
N^^) = 202, Net'' = 554, and Nc? = 1177, in very good agreement with the values 
obtained directly from Fig. 11.51 

The curve that specifies iVcV in the {Z, N) plane defines the border of stability 
for spontaneous electron decay. In nuclei, such borders of stability against sponta- 
neous proton or neutron emission are known as nucleon drip lines [79] . For the case 
of sodium clusters, the electron drip line is displayed in the inset of Fig. 11.51 
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Fig. 1.6. Calculated smooth total energy per atom as a function of the excess negative charge 
Z for the three families of sodium clusters with N = 30, N = 80, and N = 240 atoms. A 
spherical jellium background was used. As the straight lines in the inset demonstrate, the curves 
are parabolic. 'We find that they can be fitted by Eq. 111.2911 . See text for an explanation of how 
the function g{N, Z) was extracted from the calculations. 

1.3.1.3. Critical sizes for potassium and aluminum 

While in this investigation -we have used sodium clusters as a test system, the 
methodology and conclusions extend to other materials as 'well. Thus given a cal- 
culated or measured bulk 'work function W , and a spill-out parameter {5q typically 
of the order of 1-2 a.u., and neglecting 82), one can use Eq. (|1.28p . 'with Az = 0, to 
predict critical sizes for other materials. For example, our calculations for potassium 
(rs = 4.86 a.u.) give fitted values W = 2.6 eV (compared to a KS-DFT value of 
2.54 eV for a semi-infinite planar surface 'with r^. = 5.0 a.u. [78]) and 5q — 1.51 a.u. 
for (52 = 0, yielding N^V = 33, TVir^ = 152, and ivir^ =421. 

As a further example, ■we give our results for a trivalent metal, i.e. aluminum 
(rs = 2.07 a.u.), for -which our fitted values are W = 3.65 eV (compared to a 
KS-DFT value of 3.78 eV for a semi-infinite plane surface, -with = 2.0 a.u. [78]) 
and ,5o = 1.86 a.u. for 82 = 0, yielding Tvi^^ = 40 (121 electrons), ivir^ = 208 (626 
electrons), and TVir^ = 599 (1796 electrons). 

1.3.1.4. Metastability against electron autodetachment 

The multiply charged anions -with negative affinities do not necessarily exhibit a 
positive total energy. To illustrate this point, we display in Fig. 11.61 the calculated 
total energies per atom {E{N, Z)/N) as a function of excess charge {Z) for clus- 
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ters containing 30, 80, and 240 sodium atoms. These sizes allow for exothermic 
attachment of maximum one, two, or three excess electrons, respectively. 

As was the case with the electron affinities, the total-energy curves in Fig. 11.61 
show a remarkable regularity, suggesting a parabolic dependence on the excess 
charge. To test this conjecture, we have extracted from the calculated total ener- 
gies the quantity g{N, Z) = G{N, Z)/N where G(7V, Z) = [E{N, Z) - E{N, 0)]/Z + 
Ai{N), and have plotted it in the inset of Fig. 11.61 as a function of the excess neg- 
ative charge Z. The dependence is linear to a remarkable extent; for Z — 1 all 
three lines cross the energy axis at zero. Combined with the results on the electron 
affinities, this indicates that the total energies have the following dependence on 
the excess number of electrons [Z): 



where the dependence on the number of atoms in the cluster is not explicitly indi- 
cated. 

This result is remarkable in its analogy with the classical image-charge result 
of van Staveren et al. [77]. Indeed, the only difference amounts to the spill-out 
parameter 5q and to the weak dependence on Z through 62- This additional Z- 
dependence becomes negligible already for the case of 30 sodium atoms. 

For metastable multiply-charged cluster anions, electron emission (autodetach- 
ment) will occur via tunneling through a barrier (shown in Fig. ll.7p . However, to re- 
liably estimate the electron emission, it is necessary to correct the LDA effective po- 
tential for self-interaction effects. We performed a self-interaction correction of the 
Amaldi type [73] for the Hartree term and extended it to the exchange-correlation 
contribution to the total energy as follows: E^^[p] = E]^^^[p] - NeE]^^^[p/Ne], 
where N,, ~ vN + Z is the total number of electrons. This self-interaction correction 
is akin to the orbitally-averaged-potential method [73] . Minimizing the SIC energy 
functional for the parameters tq, a, and 7, we obtained the effective SIC poten- 
tial for Na^^ shown in Fig. II. 7[ which exhibits the physically correct asymptotic 
behavior [80]. 

The spontaneous electron emission through the Coulombic barrier is analogous 
to that occurring in proton radioactivity from neutron-deficient nuclei [81], as well 
as in alpha-particle decay. The transition rate is A = ln2/rx/2 = i^P, where v is 
the attempt frequency and P is the transmission coefficient calculated in the WKB 
method (for details, cf. Ref. [81]). For the 2s electron in Na^^ (cf. Fig. II. 7|) . we 
find V = 0.73 10^5 Hz and P = 4.36 10"^, yielding T1/2 = 2.18 10"!° s. For a 
cluster size closer to the drip line (see Fig. II. 5p . e.g. Nag^, we find T1/2 = 1-13 s. 

Finally, the exression in Eq. ()1.29|) for the total energy can be naturally extended 
to the case of multiply positively charged metal clusters by setting Z = —z, with 
z > 0. The ensuing equation retains the same dependence on the excess positive 
charge z, but with the negative value of the first affinity, —Ai, replaced by the 
positive value of the first ionization potential, Ii — W + (3/8)e^/(i? + S), a result 



E{Z) ^ E{0) - AiZ + 



Z{Z-l)e^ 
2{R + S) 



(1.29) 
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Fig. 1.7. The DFT (LDA) and the corresponding self-interaction corrected potential for the 
metastable Na^g" cluster. A spherical jellium background was used. The single-particle levels 
of the SIC potential are also shown. Unlike the LDA, this latter potential exhibits the correct 
asymptotic behavior. The 2s and Id electrons can be emitted spontaneously by tunneling through 
the Coulombic barrier of the SIC potential. Distances in units of the Bohr radius, ao- The specified 
single-particle levels are associated with the SIC potential. 

that has been suggested from earher measurements on multiply charged potassium 
cations [82]. Naturally, the spill-out parameter 5 assumes different values than in 
the case of the anionic clusters. 

1.3.2. Neutral and multiply charged fullerenes 

1.3.2.1. Stabilized jellium approximation - The generalized DFT-SC'M 

Fullerenes and related carbon structures have been extensively investigated using ab 
initio density-functional-theory methods and self-consistent solutions of the Kohn- 
Sham (KS) equations [83, 84]. For metal clusters, replacing the ionic cores with a 
uniform jellium background was found to describe well their properties within the 
KS-DFT method [58]. Motivated by these results, several attempts to apply the 
jellium model in conjunction with DFT to investigations of fullerenes have appeared 
recently [39, 85-87]. Our approach [39] differs from the earlier ones in several aspects 
and, in particular, in the adaptation to the case of finite systems of the stabilized- 
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jellium (or structureless pseudopotential) energy density functional (see Eq. (|1.30p 
below and Ref. [73]). 

An important shortcoming of the standard jellium approximation for fuUerenes 
(and other systems ■with high density, i.e., small Ts) results from a "well- known prop- 
erty of the jellium at high electronic densities, namely that the jellium is unstable 
and yields negative surface-energy contribution to the total energy [73], as well as 
unreliable values for the total energy. These inadequacies of the standard jellium 
model can be rectified by pseudopotential corrections. A modified-jellium approach 
which incorporates such pseudopotential corrections and is particularly suited for 
our purposes here, is the structureless pseudopotential model or stabilized jellium 
approximation developed in Ref. [73]. 

In the stabilized jellium, the total energy Epseudo, as a functional of the electron 
density p{r), is given by the expression 



Epseudolp, P+] = Ejeii[p,p+] + {Sv)ws / p{r)U{r)dr -e p+{r)dr, (1.30) 



where by definition the function U (r) equals unity inside, but vanishes, outside the 
jellium volume. p+ is the density of the positive jellium background (which for 
the case of Cgo is taken as a spherical shell, of a certain width 2d, centered at 
6.7 a.u. ). Epseudo in Eq. ()1.30|) is the standard jellium-model total energy, Ejeii, 
modified by two corrections. The first correction adds the effect of an average (i.e., 
averaged over the volume of a Wigner-Seitz cell) difference potential, {6v)wsl^{i^)^ 
which acts on the electrons in addition to the standard jellium attraction and is 
due to the atomic pseudopotentials (in this work, we use the Ashcroft empty-core 
pseudopotential, specified by a core radius rc, as in Ref. [73]). The second correction 
subtracts from the jellium energy functional the spurious electrostatic self-repulsion 
of the positive background within each cell; this term makes no contribution to the 
effective electronic potential. 

Following Ref. [73], the bulk stability condition (Eq. (25) in Ref. [73]) determines 
the value of the pseudopotential core radius Tc, as a function of the bulk Wigner- 
Seitz radius r^. Consequently, the difference potential can be expressed solely as a 
function of as follows (energies in Ry, distances in a.u.): 



where £c is the per particle electron-gas correlation energy (in our calculation, we 
use the Gunnarsson-Lundqvist exchange and correlation energy functionals; see 
Refs. [3, 23]). 

The electrostatic self-energy, e, per unit charge of the uniform positive jellium 
is given by 




(1.31) 



S I 



(1.32) 



where v is the valence of the atoms (w = 4 for carbon) . 
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1.3.2.2. ETF electron- density profile 

To apply the ETF-DFT method to carbon fullerenes, we generahze it by employing 
potential terms according to the stabilized-jellium functional in Eq. (|1.30p . 

Another required generalization consists in employing a parametrized electron- 
density profile that accounts for the hollo"w cage-structure of the fullerenes. Such 
a density profile is provided by the follo'wing adaptation of a generalization of an 
inverse Thomas-Fermi distribution, used earlier in the context of nuclear physics 
[88], i.e., 

\cosh[wi^o/ai^o\ + cosh[(r - R)/ai^o\J 

■where R = 6.7 a.u. is the radius of the fullerene cage, w, a, and 7 are variables to 
be determined by the ETF-DFT minimization. For R = and large values of w/a, 
expression (ll.33|) approaches the more familiar inverse Thomas-Fermi distribution, 
•with w the ■width, a the diffuseness and 7 the asymmetry of the profile around 
r — w. There are a total of six parameters to be determined, since the indices 
(i,o) stand for the regions inside (r < R) and outside (r > R) the fullerene cage. 
Fi^o = {cosh[w'Lo/ai.o] + 1)/ sinh[wi, o/a-i.o] is a constant guaranteeing that the t^wo 
parts of the curve join smoothly at r = i?. The density profile in Eq. (|1.33p peaks 
at r = i? and then falls to^wards smaller values both inside and outside the cage (see 
top panel of Fig. 11.81) . 

1.3.2.3. Shell correction and icosahedral splitting 

To apply the SCM to the present case, the potential Vetf in Eq- (|1.19P is replaced 
by the stabilized-jellium LDA potential sho^wn in Fig. 11.81 After some rearrange- 
ments, the shell-corrected total energy Esh[p\ in the stabilized-jellium case can be 
■written in functional form as foUo^ws [compare to Eq. (|1.2ip . see also Eq. (|1.16p ]. 

£^cW)]dr + Ej-eJ p+(r)dr, (1.34) 

Heretofore, the point-group icosahedral symmetry of Cgo ■was not considered, 
since the molecule ■was treated as a spherically symmetric cage. This is a reasonable 
zeroth-order approximation as noticed by several authors [83, 87, 89, 90]. Ho^wever, 
considerable improvement is achieved ■when the effects of the point-group icosahedral 
symmetry are considered as a next-order correction (mainly the lifting of the angular 
momentum degeneracies) . 

The method of introducing the icosahedral splittings is that of the crystal field 
theory [91]. Thus, we ■will use the fact that the bare electrostatic potential from 
the ionic cores, considered as point charges, acting upon an electron, obeys the 
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Fig. 1.8. Bottom panel: The stabilized-jellium LDA potential obtained by the ETF method 
for the neutral Ceo molecule. The Wigner-Seitz radius for the jellium bacground is 1.23 a.u. 
Note the asymmetry of the potential about the minimum. The associated difference potential 
{Sv)ws = -9.61 eV. 

Top panel: Solid line: Radial density of the positive jeUium background. Dashed Hne: ETF 
electronic density. Note its asymmetry about the maximum. Thick solid line: The difference 
(multiplied by 10) of electronic ETF densities between Cg^ and Ceo- It illustrates that the 
excess charge accumulates in the outer perimeter of the total electronic density. All densities are 
normalized to the density of the positive jellium background. 



well-known expansion theorem [91] 

U{r) = -ve-'J2T^^=-T. E ^iir)CrYr{0A), (1-35) 



oo I 



1=0 rn=-l 

where the angular coefficients C™ are given through the angular coordinates 9i,i 
of the carbon atomic cores, namely, 



cr = Y.Yr{ei,ci>i), 

i 

where * denotes complex conjugation. 



(1.36) 
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Fig. 1.9. (a) The single-particle levels of the ETF-LDA potential for Cgo shown in Fig. 11.81 
Because of the spherical symmetry, they are characterized by the two principle quantum numbers 
Ur and I, where Ur is the number of radial nodes and I the angular momentum. They are grouped 
in three bands labeled cr (rir = 0), tt {ur = 1), and S (ur = 2). Each band starts with an Z = 
level. 

(b) The single-particle levels for Ceo after the icosahedral splittings are added to the spectra of (a). 
The tenfold degenerate HOMO (h„) and the sixfold degenerate LUMOl (ti„) and LUM02 (tig) 
are denoted; they originate from the spherical 1 = 5 and 1 = 6 (tig) tt levels displayed in panel 
(a). For the a electrons, the icosahedral perturbation strongly splits the I = 9 level of panel (a). 
There result five sublevels which straddle the cr-electron gap as follows: two of them (the eightfold 
degenerate gu, and the tenfold degenerate h^) move down and are fully occupied resulting in a 
shell closure (180 cr electrons in total). The remaining unoccupied levels, originating from the 
I = 9 a level, are sharply shifted upwards and acquire positive values. 



We take the radial parameters K;(r) as constants, and determine their value by 
adjusting the icosahedral single-particle spectra e]'^" to reproduce the pseudopoten- 
tial calculation of Ref. [83] , "which are in good agreement "with experimental data. 
Our spectra ■without and "with icosahedral splitting are shown in Fig. ILOf a) and 
Fig. ll.Qr b). respectively. 

We find that a close reproduction of the results of ab initio DFT calcula- 
tions [83, 92, 93] is achieved when the Wigner-Seitz radius for the jellium back- 
ground is 1.23 a.u. The shell corrections, AE^'j", including the icosahedral splittings 
are calculated using the icosahedral single-particle energies e^"^" in Eq. ()1.19|) . The 
average quantities (p and V) are maintained as those specified through the ETF 
variation "with the spherically symmetric profile of Eq. (|1.33|) . This is because the 
first-order correction to the total energy (resulting from the icosahedral perturba- 
tion) vanishes, since the integral over the sphere of a spherical harmonic Yj"^ (l > 0) 
vanishes. 
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1.3.2.4. Ionization potentials and electron affinities 

Having specified the appropriate Wigner-Seitz radius Vs and the parameters m of 
the icosahedral crystal field through a comparison with the pseudopotential DFT 
calculations for the neutral Ceo, we calculate the total energies of the cationic and 
anionic species by allo"wing for a change in the total electronic charge, namely by 
imposing the constraint 



•where p{r) is given by Eq. (|1.33p . The shell-corrected and icosahedrally perturbed 
first and higher ionization potentials Z^'^" are defined as the difference of the ground- 
state shell-corrected total energies as follows: 



= Ell^iNe = 240 - x; Z = 240) - EH^iNe = 240 - a; + 1; Z = 240), (1.38) 



where iVe is the number of electrons in the system and x is the number of excess 
charges on the fullerenes (for the excess charge, we will find convenient to use two 
different notations x and z related as a; = |z|. A negative value of z corresponds to 
positive excess charges). Z = 240 denotes the total positive charge of the jellium 
background. 

The shell-corrected and icosahedrally perturbed first and higher electron affini- 
ties A^Jr" are similarly defined as 



^ico ^ El'j°{Ne ^2A0 + x-l;Z^ 240) - Ell°{N^ = 240 -f x; Z = 240). (1.39) 



We have also calculated the corresponding average quantities Ix and A^, which 
result from the ETF variation with spherical symmetry (that is without shell and 
icosahedral symmetry corrections) . Their definition is the same as in Eq. (|1.38p and 
Eq. (|1.39p . but with the index sh replaced by a tilde and the removal of the index 
ico. 

In our calculations of the charged fullerene molecule, the value and the icosa- 
hedral splitting parameters (k;, see Eq. (|1.35p . and discussion below it) were taken 
as those which were determined by our calculations of the neutral molecule, dis- 
cussed in the previous section. The parameters which specify the ETF electronic 
density (Eq. ()1.33|) ) are optimized for the charged molecule, thus allowing for re- 
laxation effects due to the excess charge. This procedure is motivated by results 
of previous electronic structure calculations for Cfig and Cgg [92, 93], which showed 
that the icosahedral spectrum of the neutral Ceo shifts almost rigidly upon charging 
of the molecule. 

Shell-corrected and ETF calculated values of ionization potentials and electron 
affinities, for values of the excess charge up to 12 units, are summarized in Table [LT] 
(for Vs — 1.23 a.u.) 




(1.37) 
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Table 1.1. ETF (spherically averaged, de- 
noted by a tilde) and shell-corrected (de- 
noted by a superscript ico to indicate that 
the icosahedral splittings of energy lev- 
els have been included) IPs and EAs of 
fuUerenes Cg^. Energies in eV. rs = 1.23 













X 




jico 




A ico 


1 


5.00 


7.40 


2.05 


2.75 


2 


7.98 


10.31 


-0.86 


-0.09 


3 


10.99 


13.28 


-3.75 


-2.92 


4 


14.03 


16.25 


-6.60 


-5.70 


5 


17.09 


19.22 


-9.41 


-8.41 


6 


20.18 


22.20 


-12.19 


-11.06 


7 


23.29 


25.24 


-14.94 


-14.85 


8 


26.42 


28.31 


-17.64 


-17.24 


9 


29.57 


31.30 


-20.31 


-19.49 


10 


32.73 


34.39 


-22.94 


-21.39 


11 


35.92 


39.36 


-25.53 


-22.93 


12 


39.12 


42.51 


-28.07 


-23.85 



1.3.2.5. Charging energies and capacitance of fuUerenes 

Figure ll.lOf a) sho'ws that the variation of the total ETF-DFT energy difference 
(appearance energies) AE{z) — E{z) — i?(0), as a function of excess charge z (|z| = 
x), exhibits a parabohc behavior. The inset in Fig. ll.lOT al exhibiting the quantity 

j,., = S£)_M + i„ ,1,0) 

plotted versus z (open squares), shows a straight line "which crosses the zero energy 
line at z = 1. As a result the total ETF-DFT energy has the form, 

E{z) = m+'-^^^^~Mz. (1.41) 

Equation (|1.4ip indicates that fuUerenes behave on the average like a capacitor 
having a capacitance C (the second term on the rhs of Eq. (11.411) corresponds to 
the charging energy of a classical capacitor, corrected for the self-interaction of the 
excess charge [3, 23]). We remark that regarding the system as a classical conductor, 
•where the excess charge accumulates on the outer surface, yields a value of C = 8.32 
a.u. (that is the outer radius of the jellium shell). Naturally, the ETF calculated 
value for C is somewhat larger because of the quantal spill-out of the electronic 
charge density. Indeed, from the slope of g(z) we determine [94] C — 8.84 a.u. 

A similar plot of the shell-corrected and icosahedrally modified energy differences 
AEil°{z) = Eil°{z)-Eil°{0) is shown in Fig.lHOlJb) (in the range -2 < z < 4, filled 
circles). The function gl'j^{z), defined as in Eq. (|1.40p but with the shell-corrected 
quantities (Ai?|.™(z) and Al'^"), is included in the inset to Fig. ll.lOf a) (filled circles). 
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Fig. 1.10. (a) ETF-DFT total energy differences (appearance energies) Ai?(z) = E(^z) — i?{0) as 
a function of the excess charge 2 < corresponds to positive excess charge). Inset: The ETF 
function 9(2:) (open squares), and the shell-corrected function g^^^(z) (filled circles). For z>_ \ the 
two functions are almost identical. 

(b) magnification of the appearance-energy curves for the region —2 < z < 4. Filled circles: shell- 
corrected icosahedral values \C\El^°{z) = E^^°{z) - S»^°(0)]. Open squares: ETF-DFT values 
{l^E{z) = E{z)-Em. 



The shift discernible between 1) and g^^^(X) approximately 1.7 eF, and 

originates from the difference of shell effects on the IPs and EAs (see Table [O]) . The 
segments of the curve gl'j°{z) in the inset of Fig. fTTTUT a) . corresponding to positively 
[z < 0) and negatively (z > 0) charged states, are again well approximated by 
straight lines, whose slope is close to that found for 5(2:). Consequently, we may 
approximate the charging energy, including shell-effects, as follows. 



i;r(x) = i?r(o) + 

for negatively charged states, and 



x{x — l)e^ 
2C 



x{x — l)e^ 
2C 



A\^°x, 



irx, 



(1.42) 



(1.43) 
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for positively charged states. Note that without shell-corrections (i.e., ETF only) 
7i - li = eVC = 27.2/8.84 eV « 3.1 eV, because of the symmetry of Eq. ([LiT]) 
■with respect to z, while the shell-corrected quantities are related as /J"^" — A^^" « 
e^/C -\- Ash, ■where the shell correction is Ash ~ 1-55 eV (from Table fO] 11'^° — 
° « 4.65 eV). 

Expression (|1.42|) for the negatively charged states can be rearranged as foUo^ws 
(energies in units of eV), 

Ell%x) - Ell°iO) - -2.99 + 1.54(a; - 1.39)2, (1.44) 

in close agreement ■with the all-electron LDA result of Rcf. [95] . 

Equations (fTiS]) and (fOSl) can be used to provide simple analytical approxima- 
tions for the higher IPs and EAs. Explicitly written, A^"^" = Eif{x - 1) - Ef^{x) = 
Af" ~ {x ~ l)eVC and 1^^° = /j™ + {x - l)e^/C. Such expressions have been 
used previously [96] with an assumed value for C ~ 6.7 a.u. (i.e., the radius of 
the Ceo molecule, as determined by the distance of carbon nuclei from the cen- 
ter of the molecule), which is appreciably smaller than the value obtained by us 
(C = 8.84 a.u., see above) via a microscopic calculation. Consequently, using the 
above expression with our calculated value for Af° = 2.75 eV (see Table fOT) . we 
obtain an approximate value of A^2° = —0.35 eV (compared to the microscopi- 
cally calculated value of —0.09 eV given in Table [TTTl and —0.11 eV obtained by 
Ref. [95]) — indicating metastability of CgQ — while employing an experimental 
value for A^^"" = 2.74 eV , a value of A^° = 0.68 eV was calculated in Ref. [96]. 

Concerning the cations, our expression (|1.43p with a calculated P^" = 7.40 eV 
(see Table fTTTj) and C — 8.84 a.u. yields approximate values 18.5 eV and 31.5 eV for 
the appearance energies of CgJ and Cg([ (compared to the microscopic calculated 
values of 17.71 eV and 30.99 eV , respectively, extracted from Table [TTTl and 18.6 
eV for the former obtained in Ref. [92]). Employing an experimental value for 
jico _ y ^ corresponding values of 19.20 eV and 34.96 eV were calculated in 

Ref. [96] . As discussed in Ref. [97] , these last values are rather high, and the origin 
of the discrepancy may be traced to the small value of the capacitance which was 
used in obtaining these estimates in Ref. [96]. 

A negative value of the second affinity indicates that C^^ is unstable against elec- 
tron autodetachment. In this context, we note that the doubly negatively charged 
molecule CgQ has been observed in the gas phase and is believed to be a long-lived 
metastable species [98, 99]. Indeed, as we discuss in the next section, the small DFT 
values of A'^2° found by us and by Ref. [95] yield lifetimes which are much longer 
than those estimated by a pseudopotential-like Hartree-Fock model calculation [98] , 
where a value of ~ 1 /xs was estimated. 

1.3.2.6. Lifetimes of metastable anions, Cqq~ 

The second and higher electron affinities of Ceo "were found to be negative, which 
implies that the anions Cqq with x > 2 are not stable species, and can lower their 
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energy by emitting an electron. Ho'wever, unless the number of excess electrons 
is large enough, the emission of an excess electron involves tunneling through a 
barrier. Consequently, the moderately charged anionic fullerenes can be described 
as metastable species possessing a decay lifetime. 

To calculate the lifetime for electron autodetachmant, it is necessary to deter- 
mine the proper potential that the emitted electron sees as it leaves the molecule. 
The process is analogous to alpha-particle radioactivity of atomic nuclei. The emit- 
ted electron 'will have a final kinetic energy equal to the negative of the correspond- 
ing higher EA. We estimate the lifetime of the decay process by using the WKB 
method, in the spirit of the theory of alpha-particle radioactivity, 'which has es- 
tablished that the main factor in estimating lifetimes is the relation of the kinetic 
energy of the emitted particle to the Coulombic tail, and not the details of the 
many-body problem in the immediate vicinity of the parent nucleus. 

Essential in this approach is the determination of an appropriate single-particle 
potential that describes the transmission barrier. It is ■well kno'wn that the (DPT) 
LDA potential posseses the -wrong tail, since it allo-ws for the electron to spuriously 
interact 'with itself. A more appropriate potential 'would be one produced by the Self- 
Interaction Correction method of Ref. [78] . This potential has the correct Coulombic 
tail, but in the case of the fullerenes presents another dra'wback, namely Koopman's 
theorem is not satisfied to an extent adequate for calculating lifetimes [100]. In this 
context, 'we note that Koopman's theorem is kno'wn to be poorly satisfied for the 
case of fullerenes even in Hartree-Fock calculations [101]. Therefore, the HOMO 
corresponding to the emitted electron, calculated as described above, cannot be 
used in the WKB tunneling calculation. 

Since the final energy of the ejected electron equals the negative of the value 
of the electron affinity, we seek a potential that, together -with the icosahedral 
perturbation, yields a HOMO level in C^q 'with energy —A^^°. We construct this 
potential through a self-interaction correction to the LDA potential as foUo'ws, 



where the parameter ^ is adjusted so that the HOMO level of Cqq equals —A't^°. In 
the above expression, the second term on the rhs is an average self-interaction 
Hartree correction 'which ensures a proper long-range behavior of the potential 
(i.e., correct Coulomb tail), and the third term is a correction to the short-range 
exchange-correlation. 

For the cases of Cgg" and G^^ such potentials are plotted in Fig. 11.111 We ob- 
serve that they have the correct Coulombic tail, namely a tail corresponding to one 
electron for CgQ and to t'wo electrons for C^q . The actual barrier, ho-wever, through 
'which the electron tunnels is the sum of the Coulombic barrier plus the contribution 
of the centrifugal barrier. As seen from Fig. 11.111 the latter is significant, since the 
HOMO in the fullerenes possesses a rather high angular momentum {I — 5), -while 
being confined in a small volume. 



VwKB = Vlda[p\ - Vh[4-] - V.ci^-^l 



(1.45) 
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Fig. 1.11. WKB eflfective barriers used to estimate lifetimes for Cqq (a) and Cgg" (b). Dashed 
lines correspond to barriers due solely to Coulombic repulsion and solid lines to total barriers after 
adding the centrifugal components. The thick horizontal solid lines correspond to the negative 
of the associated electron aiRnities Aj^" (a) and (b). In the case of Cg^ [panel (a)], the 

horizontal solid line at —A^^" = 0.09 eV crosses the total barrier at an inside point i?i = 9.3 a.u. 
and again at a distance very far from the center of the fuUcrene molecule, namely at an outer point 



R2 



27.2/0.09 a.u. = 302.2 a.u. This large value of R2, combined with the large 



centrifugal barrier, yields a ma<;roscopic lifetime for the metastable Cgg (see text for details). 



Using the WKB approximation [102], we estimate for C^^ a macroscopic half-life 
of ~ 4 X 10'' years, while for Cgg we estimate a very short half-life of 2.4 x 10~^^ s. 
Both these estimates are in correspondence with observations. Indeed, Cgg has 
not been observed as a free molecule, while the free Cg,^ has been observed to be 
long lived [98, 99] and was detected even 5 min after its production through laser 
vaporization [99]. 

We note that the WKB lifetimes calculated for tunneling through Coulombic 
barriers are very sensitive to the final energy of the emitted particle and can vary 
by many orders of magnitude as a result of small changes in this energy, a feature 
well known from the alpha radioctivity of nuclei [102]. 

Since the second electron affinity of Ceo is small, effects due to geometrical relax- 
ation and spin polarization can influence its value and, consequently, the estimated 
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lifetime. Nevertheless, as shown in Ref. [95], inclusion of such corrections yields 
again a negative second affinity, but of somewhat smaller magnitude, resulting in 
an even longer lifetime (the sign conventions in Ref. [95] are the opposite of ours). 

Furthermore, as discussed in Ref. [103], the stabilization effect of the Jahn- Teller 
relaxation for the singly- charged ion is only of the order of 0.03 - 0.05 eV. Since 
this effect is expected to be largest for singly-charged species, Cgp is not expected 
to be influenced by it [95]. 

On the other hand, generalized exchange-correlation functionals with gradient 
corrections yield slightly larger values for the second electron affinity. For example, 
using exchange-correlation gradient corrections, Ref. [95] found A'2° = —0.3 eV, 
which is higher (in absolute magnitude) than the value obtained without such cor- 
rections. This value of —0.3 eV leads to a much smaller lifetime than the several 
million of years that correspond to the value of —0.09 eV calculated by us. Indeed, 
using the barrier displayed in Fig. ll.llT a). we estimate a lifetime for Cgo" of approx. 
0.37 s, when Ai^" — —0.3 eV. We stress, however, that even this lower-limit value 
still corresponds to macroscopic times and is 5 orders of magnitude larger than the 
estimate of Ref. [98], which found a lifetime of 1 fis for Aif " = —0.3 eV, since it 
omitted the large centrifugal barrier. Indeed, when we omit the centrifugal barrier, 
we find a lifetime estimate of 1.4 /is, when A'^2° — —0.3 eV. 

1.3.3. On mesoscopic forces and quantized conductance in model 
metallic nanowires 

1.3.3.1. Background and motivation 

In this section, we show that certain aspects of the mechanical response (i.e., elon- 
gation force) and electronic transport (e.g., quantized conductance) in metallic 
nanowires can be analyzed using the DFT shell correction method, developed and 
applied previously in studies of metal clusters (see Sec. ll.2l and Sec. II. 3. l]) . Specifi- 
cally, we show that in a jellium-modelled, volume-conserving nanowire, variations of 
the total energy (particularly terms associated with electronic subband corrections) 
upon elongation of the wire lead to self- selection of a sequence of stable "magic" 
wire configurations (MWC's, specified in our model by a sequence of the wire's 
radii), with the force required to elongate the wire from one configuration to the 
next exhibiting an oscillatory behavior. Moreover, we show that due to the quan- 
tized nature of electronic states in such wires, the electronic conductance varies in 
a quantized step- wise manner (in units of the conductance quantum go — 2e^/h), 
correlated with the transitions between MWC's and the above-mentioned force os- 
cillations. 

Prior to introducing the model, it is appropriate to briefly review certain pre- 
vious theoretical and experimental investigations, which form the background and 
motivation for this study of nanowires. Atomistic descriptions, based on realistic 
interatomic interactions, and/or first-principles modelling and simulations played 
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an essential role in discovering the formation of nanowires, and in predicting and 
elucidating the microscopic mechanisms underlying their mechanical, spectral, elec- 
tronic and transport properties. 

Formation and mechanical properties of intcrfacial junctions (in the form of crys- 
talline nanowires) have been predicted through early molecular-dynamics simula- 
tions [104], where the materials (gold) were modelled using semiempirical embedded- 
atom potentials. In these studies it has been shown that separation of the contac;t 
between materials leads to generation of a connective junction which elongates and 
narrows through a sequence of structural instabilities; at the early stages, elonga- 
tion of the junction involves multiple slip events, while at the later stages, when the 
lateral dimension of the wire necks down to a diameter of about 15 A, further elon- 
gation involves a succession of stress accumulation and fast relief stages associated 
with a sequence of order-disorder structural transformations localized to the neck 
region [104-106]. These structural evolution patterns have been shown through the 
simulations to be portrayed in oscillations of the force required to elongate the wire, 
with a period approximately equal to the interlayer spacing. In addition, the "saw- 
toothed" character of the predicted force oscillations (see Fig. 3(b) in Ref. [104] and 
Fig. 3 in Ref. [105]) reflects the stress accumulation and relief stages of the elonga- 
tion mechanism. Moreover, the critical resolved yield stress of gold nanowires has 
been predicted [104, 105] to be ~ 4GPa, which is over an order of magnitude larger 
than that of the bulk, and is comparable to the theoretical value for Au (1.5 GPa) 
in the absence of dislocations. 

These predictions, as well as anticipated electronic conductance properties [104, 
107], have been corroborated in a number of experiments using scanning tunnel- 
ing and force microscopy [104, 108-113], break junctions [114], and pin-plate tech- 
niques [105, 115] at ambient environments, as well as under ultrahigh vacuum and/or 
cryogenic conditions. Particularly, pertinent to this section are experimental ob- 
servations of the oscillatory behavior of the elongation forces and the correlations 
between the changes in the conductance and the force oscillations; see especially the 
simultaneous measurements of force and conductance in gold nanowires in Ref. [112], 
where in addition the predicted "ideal" value of the critical yield stress has also been 
measured (see also Ref. [113]). 

The jellium-based model introduced in this paper, which by construction is 
devoid of atomic crystallographic structure, does not address issues pertaining to 
nanowire formation methods, atomistic configurations, and mechnanical response 
modes (e.g., plastic deformation mechanisms, interplanar slip, ordering and dis- 
ordering mechanisms (see detailed descriptions in Refs. [104, 105] and [106], and 
a discussion of conductance dips in Ref. [110]), defects, mechanichal reversibil- 
ity [105, 112]. and roughening of the wires' morphology during elongation [106]), 
nor does it consider the effects of the above on the electron spectrum, transport 
properties, and dynamics [116]. Nevertheless, as shown below, the model offers 
a useful framework for linking investigations of solid-state structures of reduced 
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dimensions (e.g., nano"wires) "with methodologies developed in cluster physics, as 
■well as highlighting certain nanowire phenomena of mesoscopic origins and their 
analogies to clusters. 

1.3.3.2. The jellium model for metallic nanowires: Theoretical method and 
results 

Consider a cylindrical jellium "wire of length L, having a positive background "with 
a circular cross section of constant radius i? ^ L [117]. For simplicity, we restrict 
ourselves here to this symmetry of the "wire cross section. Variations in the shape of 
the nanowire cross section serve to affect the degeneracies of the electronic spectrum 
[118, 119] "without affecting our general conclusions. We also do not include here 
variations of the "wire's shape along its axis. Adiabatic variation of the -wire's axial 
shape introduces a certain amount of smearing of the conductance steps through 
tunnelling, depending on the axial radius of curvature of the -wire [118-120]. Both 
the cross-sectional and axial shape of the "wire can be included in our model in a 
rather straightfor"ward manner. 

As elaborated in Sec. 11.21 the principal idea of the SCM is the separation of the 
total DFT energy Et{R) of the nano"wire as 

Et{R) = E{R) + AEsh{R), (1.46) 

■where E{R) varies smoothly as a function of the radius R of the "wire (instead of 
the number of electrons N used in Sec. 11.21) . and AEsh{R) is the shell-correction 
term arising from the discrete quantized nature of the electronic levels. Again, as 
elaborated in Sec. II. 2[ the smooth contribution in Eq. (jl.46l) is identified ■with 
Eetf[p\- The trial radial lateral density p{r) is given by Eq. (|1.25p . and the con- 
stant po at a given radius R is obtained under the normalization condition (charge 
neutrality) 27r / p{r)rdr = p^^\r), ■where p^^\r) — 3i?^/(4rf) is the linear posi- 
tive background density. Using the optimized p, one solves for the eigenvalues of 
the Hamiltonian H = — (?i^/2to)V^ -I- Vetf[p\, and the shell correction is given by 

AEsh = -^Harris [p\ — EeTF [p] 
occ ^ 

j p{'r)VETF[p{v)]dY - Tetf[p], (1.47) 

■where the summation extends over occupied levels. Here the dependence of all 
quantities on the pertinent size variable (i.e., the radius of the ■wire R) is not sho^wn 
explicitly. Additionally, the index i can be both discrete and continuous, and in the 
latter case the summation is replaced by an integral (see belo^w). 

Follo^wing the above procedure ■with a uniform background density of sodium 
{vs — 4 a.u.), a typical potential VETF{r) for R — 12.7 a.u., ■where r is the radial 
coordinate in the transverse plane, is sho^wn in Fig. 11.121 along ■with the transverse 
eigenvalues and the Fermi level; to simplify the calculations of the electronic 
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Fig. 1.12. Lower panel: The Vetf{'>') potential for a sodium wire with a uniform jellium back- 
ground of radius R = 12.7 a.u., plotted versus the transverse radial distance from the center of the 
wire, along with the locations of the bottoms of the subbands (namely the transverse eigenvalues 
f-nm', ri is the number of nodes in the radial direction plus one, and m is the azimuthal quantum 
number of the angular momentum). The Fermi level is denoted by a dashed line. Top panel: The 
jellium background volume density (dashed line) and the electronic volume density p(r) (solid line, 
exhibiting a characteristic spillout) normalized to bulk values are shown. 



spectrum, 'we have assumed (as noted above) R ^ L, "which allows us to express 
the subband electronic spectrum as 



i(fc^;i?) 



2m 



(1.48) 



•where is the electron "wave number along the axis of the "wire (z). 

As indicated earlier, taking the wire to be charge neutral, the electronic linear 
density, \R): must equal the linear positive background density, p^j^\R). The 
chemical potential (at T = the Fermi energy ep) for a wire of radius R is deter- 
mined by setting the expression for the electronic linear density derived from the 
subband spectra equal to p^j^\r), i.e.. 



2 -^^^ 



2m . , , _ 



(1.49) 



where the factor of 2 on the left is due to the spin degeneracy. The summand defines 
the Fermi wave vector for each subband, fcF,nm- The resulting variation of epiR) 
versus R is displayed in Fig. I1.13f a) , showing cusps for values of the radius where a 
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Fig. 1.13. Variation of the Fermi energy ep [shown in (a)] and of the conductance G (shown 
in (b) in units of go = 2e^/h), plotted versus the radius R, for a sodium nanowire. Note the 
coincidence of the cusps in ep with the step-rises of the conductance. The heights of the steps in 
G reflect the subband degeneracies due to the circular shape of the wire's cross section. 

new subband drops belo'w the Fermi level as R increases (or conversely as a subband 
moves above the Fermi level as R decreases upon elongation of the -wire). Using 
the Landauer expression for the conductance G in the limit of no mode mixing 
and assuming unit transmission coefficients, G{R) — goJ2n7n^i^pi^) ~'^nm{R)], 
•where Q is the Heaviside step function. The conductance of the nano"wire, sho"wn 
in Fig. II. 13T b). exhibits quantized step-wise behavior, with the step-rises coinciding 
with the locations of the cusps in eF{R), and the height sequence of the steps is 
1^0) 2^0) 2^0) IffO: reflecting the circular symmetry of the cylindrical wires' 
cross sections [107], as observed for sodium nanowires [114]. Solving for ep{R) [see 
Eq. (|1.49[) ]. the expression for the sum on the right-hand-side of Eq. (|1.47p can be 
written as 



which allows one to evaluate AEsh [Eq. ()1.47p ] for each wire radius R. Since the 
expression in Eq. (|1.50p gives the energy per unit length, we also calculate Eetf-, 
Tetf, and the volume integral in the second line of Eq. ()1.47p for cylindrical volumes 
of unit height. To convert to energies per unit volume [denoted as eriR), s{R), and 
Aesh{R)] all energies are further divided by the wire's cross-sectional area, ttR"^. 
The smooth contribution and the shell correction to the wire's energy are shown 
respectively in Fig. I1.14r a') and Fig. I1.14r b') . The smooth contribution decreases 
slowly towards the bulk value (—2.25 eV per atom [3]). On the other hand, the 





(1.50) 
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Fig. 1.14. (a-c); The smooth (a) and shell-correction (b) contributions to the total energy (c) 
per unit volume of the jellium-modelled sodium nanowire (in units of u = 10~^ eV/a.u.^), plotted 
versus the radius of the wire (in a.u.). Note the smaller magnitude of the shell corrections relative to 
the smooth contribution, (d-e): The smooth contribution (d) to the total force and the total force 
(e), plotted in units of nN versus the wire's radius. In (e), the zeroes of the force to the left of the 
force maxima occur at radii corresponding to the local minima of the energy of the wire (c). In (f), 
we reproduce the conductance of the wire (in units of go = 2e^ /h), plotted versus R. Interestingly, 
calculations of the conductance for the MWC's (i.e., the wire radii corresponding to the locations 
of the step-rises) through the Sharvin-Weyl formula, [119, 121] corrected for the finite height of 
the confining potential [121] (see lower panel of Fig. ll.12t . namely G = go (tS'/ A|, —aP/Xp) where 
S and P are the area and perimeter of the wire's cross section and Xp is the Fermi wavelength 
[Xp = 12.91 a.u. for Na) with ct = 0.1 (see Ref. [121]), yield results which approximate well the 
conductance values (i.e., the values at the bottom of the step-rises) shown in (f). 



shell corrections are much smaller in magnitude and exhibit an oscillatory behavior. 
This oscillatory behavior remains visible in the total energy [Fig. I1.14f c)] with 
the local energy minima occurring for values i?min corresponding to conductance 
plateaus. The sequence of i?min values defines the IVIWC's, that is a sequence of 
wire configurations of enhanced stability. 

From the expressions for the total energy of the wire [i.e., ileT(R), where fl = 
ttR^L is the volume of the wire] and the smooth and shell (subband) contributions 
to it, we can calculate the "elongation force" (EF), 



FriR) 



d[nsT{R)] _ 
dZ 

F{R) + AF,hiR) 



r,Jdi{R) d[AesH{R)] 
\ dL dL 



(1.51) 
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Using the volume conservation, i.e., d{TTR^L) = 0, these forces can be written 
as Ft{R) = {-KR^/2)deT{R)ldR, F{R) = {-kR^ /2)de{R)/dR, and AFshiR) = 
{ttR^ /2)d[Aesh{R)]/dR. F{R) and Ft{R) are shown in Fig. [TJM d.e). The os- 
cillations in the force resulting from the shell-correction contributions dominate. In 
all cases, the radii corresponding to zeroes of the force situated on the left of the 
force maxima coincide with the minima in the potential energy curve of the wire, 
corresponding to the MWC's. Consequently, these forces may be interpreted as 
guiding the self-evolution of the wire toward the MWC's. Also, all the local max- 
ima in the force occur at the locations of step-rises in the conductance [reproduced 
in Fig. I1.14r f')] , signifying the sequential decrease in the number of subbands below 
the Fermi level (conducting channels) as the wire narrows (i.e., as it is being elon- 
gated). Finally the magnitude of the total forces is comparable to the measured 
ones (i.e., in the nN range). 



1.4. Summary 

While it was understood rather early that the total energy of nuclei can be decom- 
posed into an oscillatory part and one that shows a slow "smooth" variation as a 
function of size, Strutinsky's seminal contribution [1] was to calculate the two parts 
from different nuclear models: the former from the nuclear shell model and the 
latter from the liquid drop model. In particular, the calculation of the oscillatory 
part (shell correction term) was enabled by employing an averaging method that 
smeared the single particle spectrum associated with a nuclear model potential. 

A semiempirical shell-correction method (referred to as SE-SCM) for metal clus- 
ters, that was developed in close analogy to the original phenomenological Struti- 
nsky approach, was presented in the Appendix, along with some applications to 
triaxial deformations and fission barriers of metal clusters. 

This chapter reviewed primarily the motivation and theory of a microscopic 
shell correction method based on density functional theory (often referred to as 
DFT-SCM and originally introduced in Ref. [3]). In developing the DFT-SCM, we 
have used for the shell correction term (arising from quantum interference effects) 
a derivation that differs from the Strutinsky methodology [1]. Instead, we have 
shown [3] that the shell correction term can be introduced through a kinetic-energy- 
type density functional [see Eq. (|1.19p and Eq. ()1.23p ]. 

The DFT-SCM is computationally advantageous, since it bypasses the self- 
consistent iteration cycle of the more familiar KS-DFT. Indeed, the DFT-SCM 
energy functional depends only on the single-particle density, and thus it belongs to 
the class of orbital-free DFT methods. Compared to previous OF-DFT approaches, 
the DFT-SCM represents an improvement in accuracy. 

Applications of the DFT-SCM to condensed-matter nanostructures, and in par- 
ticular metal clusters, fuUerenes, and nanowires, were presented in Sec. 11.31 
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Appendix A. Semi-empirical shell-correction method (SE-SCM) 

As mentioned above aheady [see, e.g., Sec. Il.l.lj . rather than proceed -with the 
microscopic route, Strutinsky proposed a method for separation of the total en- 
ergy into smooth and sheh-correction terms [see Eq. (jl.ip ] based on an averaging 
procedure. Accordingly, a smooth part, Egp, is extracted out of the sum of the 
single-particle energies [^^e Eq. (|1.6p . or equivalently Eq. ()1.16p -with p'" re- 

placed by p and by e^] by averaging them through an appropriate procedure. 
Usually, but not necessarily, one replaces the delta functions in the single-particle 
density of states by gaussians or other appropriate weighting functions. As a result, 
each single-particle level is assigned an averaging occupation number /i, and the 
smooth part Esp is formally "written as 

Esp = Y.eJ,. (A.l) 

i 

Consequently, the Strutinsky shell correction is given by 

occ 

AiJff = (A.2) 

1=1 

The Strutinsky prescription (jA.2[) has the practical advantage of using only 
the single-particle energies e^, and not the smooth density p. Taking advantage 
of this, the single-particle energies can be taken as those of an external potential 
that empirically approximates the self-consistent potential of a finite system. In the 
nuclear case, an anisotropic three-dimensional harmonic oscillator has been used 
successfully to describe the shell-corrections in deformed nuclei. 

The single-particle smooth part, Esp, ho"wever, is only one component of the 
total smooth contribution, {Ehf in the Hartree-Fock energy considered by 

Strutinsky). Indeed as can be seen from Eq. (II. 6p [or equivalently Eq. (|1.16l) ]. 

i?total«Ai?ff (A3) 

Strutinsky did not address the question of ho"w to calculate microscopically the 
smooth part E (-which necessarily entails specifying the smooth density p) . Instead 
he circumvented this question by substituting for E the empirical energies, Elum, 
of the nuclear liquid drop model, namely he suggested that 

Stotal « ASff + Eldm- (A.4) 

In applications of Eq. (jA.4|) , the single-particle energies involved in the averaging 
[see Eqs. (IA.1|) and (IA.2|) ] are commonly obtained as solutions of a Schrodinger equa- 
tion "with phenomenological one-body potentials. This last approximation has been 
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very successful in describing fission barriers and properties of strongly deformed 
nuclei using harmonic oscillator or Wood-Saxon empirical potentials. 

In the following (Sec. \A.l\ . we describe the adaptation of the SE-SCM approach 
to condensed-matter finite systems, and in particular to triaxially deformed metal 
clusters. Moreover (Sec. lA.2"|) . we will present several figures illustrating applications 
of the SE-SCM to investigations of the effects of triaxial shape-deformations on the 
properties of metal clusters [25, 33, 35-37] and to studies of large-scale deforma- 
tions and barriers in fission of charged metal clusters [26, 27, 38]. We note that the 
SE-SCM has been extended to incorporate electronic entropy effects at finite tem- 
peratures. This latter extension, referred to as finite-temperature (FT)-SE-SCM is 
not described here, but its theory can be found in Ref. [33]. 

We mention that, in addition, Strutinsky-type calculations using phenomenolog- 
ical potentials have been reported for the case of neutral sodium clusters assuming 
axial symmetry in Refs. [29-31, 122], and for the case of fission in Ref. [28]. 

A.l. Semiempirical shell-correction method for trisLxially deformed 
clusters 

A. 1.1. Liquid-drop model for neutral and charged deformed clusters 

For neutral clusters, the liquid-drop model [28, 48, 123] (LDM) expresses the smooth 
part, E, of the total energy as the sum of three contributions, namely a volume, a 
surface, and a curvature term, i.e., 

A„y" dr + cr j dS + Aa J dSn, (A.5) 

where dr is the volume element and dS is the surface differential element. The local 
curvature k is defined by the expression k — 0.5{R^]^^ + -R~i„), where Rmax and 
Rmin are the two principal radii of curvature at a local point on the surface of the 
jellium droplet which models the cluster. The corresponding coefficients can be de- 
termined by fitting the extended Thomas- Fermi (ETF)-DFT total energy Eetf[p] 
(see Sec. II.2.2[) for spherical shapes to the following parametrized expression as a 
function of the number, N, of atoms in the cluster [124], 

E^Pj^p = a^N + asN^/^ + acN^''^. (A.6) 

The following expressions relate the coefficients Ay, ct, and Ac to the correspond- 
ing coefficients, (a's), in Eq. (jA.6p . 

3 1 1 

Av = - — ^a^ ; o- = - — ttUs ; A^ = Uc- (A. 7) 

47rrj 47rr^ AttTs 

In the case of ellipsoidal shapes the areal integral and the integrated curvature 
can be expressed in closed analytical form with the help of the incomplete elliptic 
integrals J^('0, k) and £('0, k) of the first and second kind [125], respectively. Before 
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writing the formulas, we need to introduce some notations. Volume conservation 
must be employed, namely 

ab'c/Rl = abc^l, (A.8) 

where Rq is the radius of a sphere with the same volume {Rq = VgN^^^ is taken to 
be the radius of the positive jellium assuming spherical symmetry), and a = a' /Rq, 
etc., are the dimensionless semi- axes. The eccentricities are defined through the 
dimensionless semi-axes as follows 

e\ = l- [c/af 
el = l- {b/af 

el = l-{c/bf. (A.9) 

The semi-axes are chosen so that 

a>b>c. (A.IO) 

With the notation sin^ = ei, k2 = £2/^1, and fcs = es/ei, the relative (with 
respect to the spherical shape) surface and curvature energies are given [126] by 



^surf _ ab 



^sph 2 

surf 



1 - 
ei 



and 



Ecurv be 



1 + — ((1 - e?)J^(V, k2) + el£{^, fca)) 

^1 



(A.11) 



(A.12) 



The change in the smooth part of the cluster total energy due to the excess charge 
±Z was already discussed by us for spherical clusters in the previous section. The 
result may be summarized as 

AE'P''{Z) = E'P'^iZ) - E'P''{0) = tWZ -h ^^^^ ^+S)^^ ' (^'■^^^ 

where the upper and lower signs correspond to negatively and positively charged 
states, respectively, W is the work function of the metal, _Ro is the radius of the 
positive jellium assuming spherical symmetry, and 5 is a spillout-type parameter. 

To generalize the above results to an ellipsoidal shape, (t){Ro + S) = e^/ {Rq -\- 6), 
which is the value of the potential on the surface of a spherical conductor, needs 
to be replaced by the corresponding expression for the potential on the surface of 
a conducting ellipsoid. The final result, normalized to the spherical shape, is given 
by the expression 

AE''"(Z)±WZ bc^^ ^ , , 

^ ' -J'{ip,k2), (A.14) 



AE^ph{Z)±WZ ei 

where the ± sign in front of WZ corresponds to negatively and positively charged 
clusters, respectively. 
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Fig. A.l. (a) Experimental yields (denoted as "REL. ABUND.") of dianionic silver clusters Ag^y , 
plotted versus cluster size. The error bars indicate the statistical uncertainty, (b) Theoretical FT- 
SE-SCM [33] second electron affinities A2 for Ag^ clusters at T = 300 K. LDM results are 
depicted by the dashed line. The figure was reproduced from Ref. [36]. 



A. 1.2. The modified Nilsson potential 

A natural choice for an external potential to be used for calculating shell corrections 
with the Strutinsky method is an anisotropic, three-dimensional oscillator "with an 1^ 
term for lifting the harmonic oscillator degeneracies [127]. Such an oscillator model 
for approximating the total energies of metal clusters, but "without separating them 
into a smooth and a shell-correction part in the spirit of Strutinsky' s approach, has 
been used [58] "with some success for calculating relative energy surfaces and defor- 
mation shapes of metal clusters. However, this simple harmonic oscillator model 
has serious limitations, since i) the total energies are calculated by the expression 
l^jEi, and thus do not compare with the total energies obtained from the KS- 
DFT approach, ii) the model cannot be extended to the case of charged (cationic or 
anionic) clusters. Thus absolute ionization potentials, electron affinities, and fission 
energetics cannot be calculated in this model. Alternatively, in our approach, we 
are making only a limited use of the external oscillator potential in calculating a 
modified Strutinsky shell correction. Total energies are evaluated by adding this 
shell correction to the smooth LDM energies. 

In particular, a modified Nilsson Hamiltonian appropriate for metal clusters 
[128, 129] is given by 



(A.15) 
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where Hq is the hamiltonian for a three-dimensional anisotropic oscillator, namely 
Ho = A +^(a.?x2 + ujy' + ujz^) = 

Y^i^lak + -)tvjJk- (A.16) 

Uo in Eq. (|Al5l) is a dimensionless parameter, which for occupied states may 
depend on the principal quantum number n = ni+n2~\-n^ of the spherical-oscillator 
major shell associated with a given level (ni, 77,2, "-3) of the hamiltonian (for clus- 
ters comprising up to 100 valence electrons, only a weak dependence on n is found, 
see Table I in Ref. [25]). C/q vanishes for values of n higher than the corresponding 
value of the last partially (or fully) filled major shell in the spherical limit. 

1^ = X]/c=i 'fc is ^ "stretched" angular momentum which scales to the ellipsoidal 
shape and is defined as follows, 

ll = {q,P2-q2Plf, (A.17) 

(with similarly obtained expressions for li and I2 via a cyclic permutation of in- 
dices) where the stretched position and momentum coordinates are defined via the 
corresponding natural coordinates, g^'"* and as follows, 

^ qr\m,u:Uhf'^ = , (fc = 1, 2, 3), (A.18) 

= pT\l/hm,LOkf/^ = '^^^ Ak^l, 2, 3). (A.19) 

The stretched 1^ is not a properly defined angular-momentum operator, but has 
the advantageous property that it does not mix deformed states which correspond 
to sherical major shells with different principal quantum number n = ni + n2 + 
(see, the appendix in Ref. [25] for the expression of the matrix elements of P). 

The subtraction of the term < P >n— n{n + 3)/2, where < >„ denotes the 
expectation value taken over the nt/i-major shell in spherical symmetry, guaranties 
that the average separation between major oscillator shells is not affected as a result 
of the lifting of the degeneracy. 

The oscillator frequencies can be related to the principal semi-axes a', b' , and 
c' [see Eq. (IA.8I) ] via the volume-conservation constraint and the requirement that 
the surface of the cluster is an equipotential one, namely 

U!ia' = LU2b' = oj^c' = ujqRq, (A. 20) 

where the frequency ujq for the spherical shape (with radius Rq) was taken according 
to Ref. [130] to be 

(A.21) 



^ ,,,, 49eVbohr^ 

^0{N) = ^2jVl/3 



,7Vl/3 
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Fig. A. 2. Ionization potentials of neutral Kjv clusters at three temperatures, T = 10, 300, and 
500 K. Solid dots: theoretical FT-SE-SCM [33] results. Open squares: experimental measure- 
ments [129]. The best agreement between theory and experiment happens for T = 300 K (room 
temperature), indicating the importance of the electronic entropy in quenching the shell effects. 

Since in this paper we consider solely monovalent elements, N in Eq. (|A.2ip is 
the number of atoms for the family of clusters , Ts is the Wigner-Seitz radius 
expressed in atomic units, and t denotes the electronic spillout for the neutral cluster 
according to Ref. [130]. 

A. 1.3. Averaging of single-particle spectra and semi- empirical shell 
correction 

Usually Esp [see Eqs. (|A.1[) and (IA.2I) ] is calculated numerically [131]. However, 
a variation of the numerical Strutinsky averaging method consists in using the 
semiclassical partition function and in expanding it in powers of h?. With this 
method, for the case of an anisotropic, fully triaxial oscillator, one finds [43, 132] 
an analytical result, namely [133] 

X (liSN.r^' + ^^f^^^(3iVe)^/^) , (A.22) 



where iVe denotes the number of delocalized valence electrons in the cluster. 
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In the present -work, expression (jA.22[) (as modified belo'w) will be substituted 
for the average part Esp in Eq. (jA.2p . -while the sum Ei wiW be calculated 

numerically by specifying the occupied single-particle states of the modified Nilsson 
oscillator represented by the hamiltonian (jA.15[) . 

In the case of an isotropic oscillator, not only the smooth contribution, E°p'^, 
but also the Strutinsky shell correction ()A.2p can be specified analytically, [43] with 
the result 

Ai?i%(^) = ^^o(37Ve)2/3(-l + 12x(l - x)), (A.23) 

where x is the fractional filling of the highest partially filled harmonic oscillator shell. 
For a filled shell {x = 0), A£'f^*;o(0) = -^huJoi^Ne)'^^^ , instead of the essentially 
vanishing value as in the case of the ETF-DFT defined shell correction (cf. Fig. 1 of 
Ref. [25]). To adjust for this discrepancy, we add -Ai?f^*g(0) to AiJf^*'" calculated 
through Eq. (|A.2[) for the case of open-shell, as well as closed-shell clusters. 

A. 1.4. Overall procedure 

We are now in a position to summarize the calculational procedure, which consists 
of the following steps: 

(1) Parametrize results of ETF-DFT calculations for spherical neutral jellia accord- 
ing to Eq. (jA.ep . 

(2) Use above parametrization (assuming that parameters per differential element 
of volume, surface, and integrated curvature are shape independent) in Eq. 
(jA.5|) to calculate the liquid-drop energy associated with neutral clusters, and 
then add to it the charging energy according to Eq. (|A.14|) to determine the 
total LDM energy E. 

(3) Use Equations (fOs]) and (|XT6)) for a given deformation [i.e., a', b', c', or 
equivalently wi, lu2, uj3, see Eq. (jA.20p ] to solve for the single-particle spectrum 

(£.)• 

(4) Evaluate the average, Egp, of the single-particle spectrum according to Eq. 
(|A.22[) and subsequent remarks. 

(5) Use the results of steps 3 and 4 above to calculate the shell correction AiJ^^"" 
according to Eq. (|A.2p . 

(6) Finally, calculate the total energy Esh as the sum of the liquid-drop contribution 
(step 2) and the shell correction (step 5), namely Est = E + AE^il"^. 

The optimal ellipsoidal geometries for a given cluster M^^, neutral or charged, 
are determined by systematically varying the distortion (namely, the parameters a 
and b) in order to locate the global minimum of the total energy Esh {N, Z) . 
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Fig. A. 3. Two-center-oscillator [26, 27] SE-SCM results for the asymmetric channel Na^jJ" — >■ 
Na^-t-Nag". The final configuration of Na^J" is spherical. For the heavier fragment Na^", we present 
results associated with three different final shape configurations, namely, oblate [(o,s); left], spher- 
ical [(s,s); middle], and prolate [(p,s); right]. The ratio of shorter over longer axis is 0.555 for the 
oblate case and 0.75 for the prolate case. 

Bottom panel: LDM energy (surface plus Coulomb, dashed curve) and total potential energy 
(LDM plus shell corrections, solid curve) as a function of fragment separation d. The empty ver- 
tical arrow marks the scission point. The zero of energy is taken at d = 0. A number (—1.58 eV 
or —0.98 eV), or a horizontal solid arrow, denotes the corresponding dissociation energy. 
Middle panel: Shell-correction contribution (solid curve), surface contribution (upper dashed 
curve), and Coulomb contribution (lower dashed curve) to the total energy, as a function of 
fragment separation d. 

Top panel: Single-particle spectra as a function of fragment separation d. The occupied (fully or 
partially) levels are denoted with solid lines. The unoccupied levels are denoted with dashed lines. 
On top of the figure, four snapshots of the evolving cluster shapes are displayed. The solid ver- 
tical arrows mark the corresponding fragment separations. Observe that the doorway molecular 
configurations correspond to the second snapshot from the left. Notice the change in energy scale 
for the middle and bottom panels, as one passes from (o,s) to (s,s) and (p,s) final configurations. 

A. 2. Applications of SE-SCIVI to metal clusters 



As examples of applications of the SE-SCIVI, we present here three cases. In Fig. lA.Il 
we show experimental electron affinities for doubly negatively charged silver clus- 
ters [134] and compare them with theoretical calculations [36]. In Fig. lA.2i we com- 
pare FT-SE-SCIvI calculations for the IPs of neutral K^r clusters with experimental 
results [33] ; such comparisons demonstrate the importance of electronic-entropy ef- 
fects. Finally, in Fig. IA.31 we display SE-SCM calculations for the fission barriers 
associated with the asymmetric channel Na^fJ' — )■ Na^-|-Naj^ [26, 27]; see caption for 
details. The phenomenological binding potential as a function of fission-fragment 
separation is described via a two-center-oscillator model [26, 27, 135]. 

A fourth application of the SE-SCM describing the IPs of triaxially deformed 
cold sodium clusters was already used in the introductory Sec. ll.l.2l [see Fig. ll.ir c)]. 
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